---
title: "Replication R Code"
author: "Pascal Mounchid / Max Heermann"
date: "23/01/2023"
output: pdf_document
---

# 0 Preparation and data import

## 0.1 Path for the data

```{r}
setwd("")
```

## 0.2 Load packages

```{r}
library(haven)
library(tidyverse)
library(stargazer)
library(cregg)
library(scales)
library(rio)
library(Hmisc)
library(margins)
library(survey)
library(srvyr)
library(patchwork)
library(dplyr)
library(readxl)
library(openxlsx)
library(svglite)
library(sjPlot)
```

## 0.3 Import dataset

```{r}
Dat_Vig <- data.frame(as_factor(read_dta("data.dta"))) 
```


# 1 Figure 1: Support for vaccine sharing, separated by experimental groups

## 1.1 Weighting the observations and calculating relative shares

```{r}
Data_Scale <- Dat_Vig %>% 
  select(Treatment, Selection, gewichtRespondi)

Data_Scale$Selection <- as.numeric(Data_Scale$Selection) 
```

### Treatment group

```{r}
Data_Scale_Treat <- Data_Scale %>% 
  filter(Treatment == "Treatmentgroup")

mean(Data_Scale_Treat$Selection)
sd(Data_Scale_Treat$Selection)

variablen_Likert <- "Selection"
 
list_of_variables <- vector(mode = "list")
counter <- 1
 
weighted_Likert <- Data_Scale_Treat %>%
  as_survey_design(weights = gewichtRespondi)
 
for ( i in variablen_Likert) {

list_of_variables [[counter]] <- weighted_Likert %>% 
      group_by( !!ensym(i) ) %>%
      survey_tally() %>%
        rename(Variable = !!ensym(i),
              Freq = n)

counter <- counter+1
 
}

Likert_Gewicht_T <- bind_rows(list_of_variables)

Likert_Gewicht_T <- Likert_Gewicht_T %>% 
   drop_na("Variable")

Likert_Gewicht_T <- Likert_Gewicht_T %>% 
  select(!n_se) %>% 
  mutate(RelFreq = Freq * 100/sum(Freq)) 


Likert_Gewicht_T$Variable <- as.numeric(Likert_Gewicht_T$Variable)

Likert_Gewicht_T <- Likert_Gewicht_T %>% 
  add_column(Group = "Treatmentgroup")
```


### Control group

```{r}
Data_Scale_Control <- Data_Scale %>% 
  filter(Treatment == "Controlgroup")

mean(Data_Scale_Control$Selection)
sd(Data_Scale_Control$Selection)

variablen_Likert <- "Selection"
 
list_of_variables <- vector(mode = "list")
counter <- 1
 
weighted_Likert <- Data_Scale_Control %>%
  as_survey_design(weights = gewichtRespondi)
 
for ( i in variablen_Likert) {

list_of_variables [[counter]] <- weighted_Likert %>% 
      group_by( !!ensym(i) ) %>%
      survey_tally() %>%
        rename(Variable = !!ensym(i),
              Freq = n)

counter <- counter+1
 
}

Likert_Gewicht_C <- bind_rows(list_of_variables)

Likert_Gewicht_C <- Likert_Gewicht_C %>% 
   drop_na("Variable")

Likert_Gewicht_C <- Likert_Gewicht_C %>% 
  select(!n_se) %>% 
  mutate(RelFreq = Freq * 100/sum(Freq)) 


Likert_Gewicht_C$Variable <- as.numeric(Likert_Gewicht_C$Variable)

Likert_Gewicht_C <- Likert_Gewicht_C %>% 
  add_column(Group = "Controlgroup")
```


```{r}
Plot_Likert <- bind_rows(Likert_Gewicht_C, Likert_Gewicht_T)

Plot_Likert$RelFreq <- format(round(
                              Plot_Likert$RelFreq, 2), 
                              nsmall = 2)


Plot_Likert$Freq <- format(round(
                              Plot_Likert$Freq, 2), 
                              nsmall = 2)


Plot_Likert$RelFreq <- as.numeric(
                              Plot_Likert$RelFreq)


Plot_Likert$Freq <- as.numeric(
                              Plot_Likert$Freq)

Plot_Likert$Group <- factor(
                            Plot_Likert$Group,
                            levels = c('Treatmentgroup','Controlgroup'),
                                     ordered = TRUE)

Plot_Likert$Variable <- as.factor(
                                  Plot_Likert$Variable)
```

## 1.2 Figure 1

```{r}
F1 <- ggplot(Plot_Likert,                                     
                              aes(x = Variable,
                                  y = RelFreq,
                                      fill = Group)) +
                                        geom_bar(stat = "identity",
                                                 position = "dodge") +
                        scale_fill_manual(values=c("#FA8072",
                                                   "#1E90FF")) +
                                 theme(legend.position = "none") +
                                         labs(title = "Public support for vaccine donations",
                     y = "Relative share of the experimental groups", x = "Support for vaccine donations")

F1
```

```{r}

ggsave("F1.png", width = 18, height = 12, units = "cm", dpi = 300)

ggsave("F1.svg",  width = 18, height = 12, units = "cm", dpi = 300)
```

#  2 Main model: Conditional Public Support for International Vaccine Solidarity

```{r}
amces <- cj(Dat_Vig, Selection ~  Treatment + 
                                  Need + 
                                  Kontrol + 
                                  Demok + 
                                  Geo +
                                  Kosten,
            id = ~ID,
            weights = ~gewichtRespondi)

#write.xlsx(amces, 'AMCE_estiamtion.xlsx')

head(amces[c("feature", "level", "estimate", "std.error")], 20L)
```


## Figure 2

```{r}
amcesplot <- plot(amces, size = 2)

F2 <- amcesplot +
            scale_colour_grey(end = 0) + 
                            theme_bw()+
                    labs(x = "AMCE",
                        title = "Conditional support for vaccine donations") +
                                theme(text=element_text(size=17), 
                                      legend.position = 'none') +
                                      theme(plot.caption = element_text(
                                        hjust = 0, 
                                        face= "italic"), 
                                        plot.title.position = "plot", 
                                        plot.caption.position =  "plot") 
F2

```

```{r}
ggsave("F2.png", width = 19, height = 27, units = "cm", dpi = 700)

ggsave("F2.svg",  width = 18, height = 12, units = "cm", dpi = 300)
```



# 3 Linear probability models

## Dichotomize dependent variable 

```{r}
Dat_Vig <- Dat_Vig %>% 
      mutate(Dicho_DV = case_when(Selection == 1 ~ 0,
                              Selection == 2 ~ 0,
                              Selection == 3 ~ 0,
                              Selection == 4 ~ 0,
                              Selection == 5 ~ 1,
                              Selection == 6 ~ 1,
                              Selection == 7 ~ 1))


Dat_Vig$Dicho_DV <- as.numeric(Dat_Vig$Dicho_DV)
```

## Model 1

```{r}
Regress_1  <-  lm(Dicho_DV ~ Treatment + 
                             Need + 
                             Kontrol +
                             Demok + 
                             Geo + 
                             Kosten,
                             weights = gewichtRespondi, 
                             data = Dat_Vig) 

summary(Regress_1)

CI <- confint(Regress_1, level = 0.95)

CI <- as.data.frame(CI)
```

## Model 2 

### prepare socio-demographic control variables

```{r}
Dat_Vig <- Dat_Vig %>% 
 filter(!Gender == "divers")

Dat_Vig$Gender <- droplevels(Dat_Vig$Gender)


table(Dat_Vig$Gender)


Dat_Vig <- Dat_Vig %>% 
              mutate(Education = case_when(stat3_gr3 == 1 ~ "Low",
                                    stat3_gr3 == 2 ~ "Middle",
                                    stat3_gr3 == 3 ~ "High"))

Dat_Vig$Education <- as.factor(Dat_Vig$Education)

Dat_Vig$Education <- relevel(Dat_Vig$Education, ref = "Low")

table(Dat_Vig$Education)


Dat_Vig$Age <- as.numeric(Dat_Vig$Age)

table(Dat_Vig$Age)


Dat_Vig$Altruism_1 <- as.numeric(Dat_Vig$Altruism_1)
```

```{r}
Regress_2 <- lm(Dicho_DV ~ Treatment + 
                           Need + 
                           Kontrol + 
                           Demok + 
                           Geo + 
                           Kosten +
                           Age +
                           Education + 
                           Gender + 
                           Altruism_1,
                                weights = gewichtRespondi, 
                          data = Dat_Vig) 


summary(Regress_2)

CI_2 <- confint(Regress_2, level = 0.95)

CI_2 <- as.data.frame(CI_2)
```

## Supplementary Table 5  

```{r}
stargazer(Regress_1, Regress_2, type = "html", out =  "ST5.doc")
```













